!> author: 左志华
!> date: 2022-07-02
!>
!> Internal forces <br>
!> 内力：粘性力、压力
!> @todo 解决初始时刻边界压力扰动问题
module internal_force_m

    use math_m, only: r23, rp, get_distance
    use region_type_m, only: region_t
    use smoothed_kernel_function_m, only: cubic_spline_kernel
    implicit none
    private

    public :: internal_force_without_visc

contains

    !> Internal force function without viscosity <br>
    !> 压力求解
    subroutine internal_force_without_visc(region, acc, n)
        type(region_t), intent(inout) :: region !! region <br>
                                                !! 区域
        real(rp), intent(inout) :: acc(:, :)    !! accleration <br>
                                                !! 加速度
        integer, intent(in) :: n                !! ntotal <br>
                                                !! 区域内总粒子数
                                                !! @note 未实际使用
        integer :: i
        real(rp) :: h(2), w(9), dwdx(2, 9), dx(2), r

        ! --------------------------------- 计算粒子 6 的核函数与导数值 ---------------------- !
        do i = 1, 9
            call get_distance(region%loc(:, 6), region%loc(:, i), dx, r)
            call cubic_spline_kernel(r, dx, 1.0_rp, w(i), dwdx(:, i))
        end do

        ! ------------------------------------ 计算压力 ------------------------------------ !
        do i = 1, 9

            ! 用到了虚粒子的压强对实粒子的影响，防止边界处实粒子有较大的、指向边界的初始压强合力
            h(:) = -(region%pressure(6) + region%pressure(i))*dwdx(:, i)/(region%rho(6)*region%rho(i))
            acc(:, 6) = acc(:, 6) + h(:)*region%mass(i)
            print *, i, h(:)*region%mass(i)

        end do

    end subroutine internal_force_without_visc

end module internal_force_m
